poplars <- read.table('ecol 563/poplars.txt', header=T)
poplars[1:8,]
 
# split the values of trt into individual characters
out.sp <- strsplit(as.character(poplars$trt), split="")
# out.sp is a list
out.sp
out.sp[[4]]
out.sp[[47]]
out.sp[[47]][2]
 
# testing the whole list does not work
'L' %in% out.sp
 
# use sapply to test one observation at a time
sapply(out.sp, function(x) 'L' %in% x)
# alternative use of sapply
sapply(1:48, function(x) 'L' %in% out.sp[[x]])
 
# make the sapply expression the first argument of ifelse
ifelse(sapply(1:48, function(x) 'L' %in% out.sp[[x]]), 'present', 'absent')
 
# convert the result to a factor
L <- factor(ifelse(sapply(1:48, function(x) 'L' %in% out.sp[[x]]), 'present', 'absent'))
L
 
# second method using grep
# create a clean copy of poplars
poplars2 <- poplars[,1:3]
poplars2[1:4,]
# create a variable L all of whose values are 'absent'
poplars2$L <- 'absent'
poplars2[1:4,]
 
# grep is one of a family of search and replace functions
?grep
# grep returns the locations in the vector trt where 'L' is found
grep('L', poplars2$trt)
# use the locations change the value of trt to 'present'
poplars2$L[grep('L', poplars2$trt)] <- 'present'
poplars2$L
#convert the result to a factor
poplars2$L <- factor(poplars2$L)
 
# graphically test for a treatment x block interaction
interaction.plot(poplars$trt, poplars$block, poplars$weight)
 
# fit the full interaction model with blocks
out1 <- lm(weight~factor(block)+L*P*N*K, data=poplars) 
anova(out1)
 
# compare the output to the random effects version of the model
library(nlme)
lme.model <- lme(weight~L*P*N*K, random=~1|block, data=poplars) 
 
# the F-tests are nearly the same as the fixed blocks version so we can trust them
anova(lme.model)
# drop non-significant interactions
lme.model2 <- lme(weight~L*N + P + K, random=~1|block, data=poplars) 
anova(lme.model2)
# drop non-significant main effect
lme.model2 <- lme(weight~L*N + P, random=~1|block, data=poplars) 
anova(lme.model2)
fixef(lme.model2)
levels(poplars$L)
 
# create a data frame of the unique combinations of L and N
g <- expand.grid(L=levels(poplars$L), N=levels(poplars$N))
g
fixef(lme.model2)
 
# write a function to return values of the regressors for choices of L, N, and P
cvec <- function(L,N,P) c(1, L=='present', N=='present', P=='present', (L=='present')*(N=='present'))
 
# make function to obtain regressors for all combinations of L and N for a specific P
cmat <- function(P) t(apply(g, 1, function(x) cvec(x[1], x[2], P)))
 
# construct regressors when P = 'absent'
cmat1 <- cmat('absent')
# construct regressors when P = 'present'
cmat2 <- cmat('present')
cmat1
cmat2
 
# difference-adjusted confidence interval function from lecture 8
ci.func2 <- function(rowvals, glm.vmat) {
nor.func1a <- function(alpha, sig) 1-pnorm(-qnorm(1-alpha/2) * sum(sqrt(diag(sig))) / sqrt(c(1,-1) %*% sig %*%c (1,-1))) - pnorm(qnorm(1-alpha/2) * sum(sqrt(diag(sig))) / sqrt(c(1,-1) %*% sig %*% c(1,-1)), lower.tail=F)
nor.func2 <- function(a,sigma) nor.func1a(a,sigma)-.95
n <- length(rowvals)
xvec1b <- numeric(n*(n-1)/2)
vmat <- glm.vmat[rowvals,rowvals]
ind <- 1
for(i in 1:(n-1)) {
for(j in (i+1):n){
sig <- vmat[c(i,j), c(i,j)]
#solve for alpha
xvec1b[ind] <- uniroot(function(x) nor.func2(x, sig), c(.001,.999))$root
ind <- ind+1
}}
1-xvec1b
}
 
# in turth we could use a t-distribution rather than a normal distribution
names(lme.model)
# degrees of freedom are located in the $fixDF$X component
lme.model2$fixDF$X
 
# choosing t or normal makes a small difference here
qnorm(.025)
qt(.025,41)
 
# obtain variance-covariance matrices of the treatment means
vmat1 <- cmat1 %*% vcov(lme.model2) %*% t(cmat1)
vmat2 <- cmat2 %*% vcov(lme.model2) %*% t(cmat2)
 
# obtain difference-adjusted confidence intervals
ci.func2(1:4,vmat1)
ci.func2(1:4,vmat2)
 
# rewrite the make.data function from lecture 9 to work here
make.data <- function(P,model) {
#variance-covariance matrix of the means
vmat <- cmat(P) %*% vcov(model) %*% t(cmat(P))
# difference-adjusted confidence levels
vals1 <- ci.func2(1:4, as.matrix(vmat))
# levels of lh and n, means, and standard errors
part1 <- data.frame(g,
est=as.vector(cmat(P)%*%fixef(model)), se=sqrt(diag(vmat)))
# 95% intervals
part1$low95 <- part1$est + qnorm(.025)*part1$se
part1$up95 <- part1$est + qnorm(.975)*part1$se
# difference-adjusted intervals
part1$low85 <- part1$est + qnorm((1-min(vals1))/2)*part1$se
part1$up85 <- part1$est + qnorm(1-(1-min(vals1))/2)*part1$se
# add value of func and p
part1$P <- P
# return value of function
part1
}
 
# generate data for the graph
part0a <- make.data('absent',lme.model2)
part0b <- make.data('present',lme.model2)
# combine into a single data frame
fac.vals <- rbind(part0a,part0b)
fac.vals
 
# the only change needed for the group panels function is for the grid lines
 my.panel <- function(x, y, subscripts, col, pch, group.number, ...) {
#subset variables for the current panel
low95 <- fac.vals$low95[subscripts]
up95 <- fac.vals$up95[subscripts]
low85 <- fac.vals$low85[subscripts]
up85 <- fac.vals$up85[subscripts]
myjitter <- c(-.05,.05)
col2 <- c('tomato','grey60')
#95% confidence interval
panel.arrows(x+myjitter[group.number], low95, x+myjitter[group.number], up95, angle=90, code=3, length=.05, col=col)
#difference-adjusted confidence interval
panel.segments(x+myjitter[group.number], low85, x+myjitter[group.number], up85, col=col2[group.number], lineend=1, lwd=5)
# connect means with line segments
panel.lines(x+myjitter[group.number], y, col=col)
# add point estimates for the means
panel.xyplot(x+myjitter[group.number], y, col='white', pch=16)
panel.xyplot(x+myjitter[group.number], y, col=col, pch=pch, ...)
# grid lines
panel.abline(h=seq(0,70,2), col='lightgrey', lty=3)
}
 
# prepanel function needs no changes
 prepanel.ci2 <- function(x, y, ly, uy, subscripts, ...){
list(ylim=range(y, uy[subscripts], ly[subscripts], finite=T))}
 
# P is not a factor
levels(fac.vals$P)
library(lattice)
 
# modidy the xyplot call to work with the current data set
# change x-variable, conditioning variable, groups variable, x-axis labael, and the layout
xyplot(est~L|factor(P, level=c('absent','present'), labels=paste('P = ', c('absent','present'), sep='')), xlab ='L', ylab='Estimated mean', groups=N, data=fac.vals, col=c(2,1), pch=c(1,16), prepanel=prepanel.ci2, ly=fac.vals$low95, uy=fac.vals$up95, layout=c(2,1), panel.groups=my.panel, panel="panel.superpose")
 
# redo graph with a legend
xyplot(est~L|factor(P, level=c('absent','present'), labels=paste('P = ', c('absent','present'), sep='')), xlab ='L', ylab='Estimated mean', groups=N, data=fac.vals, col=c(2,1), pch=c(1,16), prepanel=prepanel.ci2, ly=fac.vals$low95, uy=fac.vals$up95, layout=c(2,1), panel.groups=my.panel, panel="panel.superpose", key=list(x=.05, y=.82, corner=c(0,0), points=list(pch=c(16,1), col=1:2, cex=.9), text=list(c('N = 1','N = 0'), cex=.8)))
 
# check to see if graphical tests coincide with summary table. The test across panels is wrong.
summary(lme.model2)
 
# generate values of L, N, and P for all 8 means together
g <- expand.grid(L=levels(poplars$L), N=levels(poplars$N), P=levels(poplars$P))
g
 
# obtain values of the regressors for each of the 8 means
cmat <- t(apply(g, 1, function(x) cvec(x[1], x[2], x[3])))
cmat
# obtain the variance-covariance matrix of the means
vmat <- cmat %*% vcov(lme.model2) %*% t(cmat)
 
# obtain difference-adjusted confidence levels
# we see they fall into three categories: .58, .75, .84
ci.func2(1:8, vmat)
 
# write version of make.data function in which the 
# difference-adjusted confidence level is an argument
make.data4 <- function(vals1,model) {
#variance-covariance matrix of the means
vmat <- cmat %*% vcov(model) %*% t(cmat)
# levels of lh and n, means, and standard errors
part1 <- data.frame(g,
est=as.vector(cmat%*%fixef(model)), se=sqrt(diag(vmat)))
# 95% intervals
part1$low95 <- part1$est + qnorm(.025)*part1$se
part1$up95 <- part1$est + qnorm(.975)*part1$se
# difference-adjusted intervals
part1$low85 <- part1$est + qnorm((1-min(vals1))/2)*part1$se
part1$up85 <- part1$est + qnorm(1-(1-min(vals1))/2)*part1$se
# return value of function
part1
}
 
# redo using .84 for the difference-adjusted confidence intervals
fac.vals <- make.data4(.84,lme.model2)
 
# generate the graph: we see that we draw the same conclusions as with .75
windows()
xyplot(est~L|factor(P, level=c('absent','present'), 
labels=paste('P = ', c('absent','present'), sep='')), 
xlab ='L', ylab='Estimated mean', groups=N, data=fac.vals, 
col=c(2,1), pch=c(1,16), prepanel=prepanel.ci2, ly=fac.vals$low95, 
uy=fac.vals$up95, layout=c(2,1), panel.groups=my.panel, panel="panel.superpose",
key=list(x=.05, y=.82, corner=c(0,0), points=list(pch=c(16,1), col=1:2, cex=.9), 
text=list(c('N = 1','N = 0'), cex=.8)))
 
# rewrite make.data function so that it returns two difference-adjusted intervals
make.data5 <- function(vals1,vals2, model) {
#variance-covariance matrix of the means
vmat <- cmat %*% vcov(model) %*% t(cmat)
# levels of lh and n, means, and standard errors
part1 <- data.frame(g,
est=as.vector(cmat%*%fixef(model)), se=sqrt(diag(vmat)))
# 95% intervals
part1$low95 <- part1$est + qnorm(.025)*part1$se
part1$up95 <- part1$est + qnorm(.975)*part1$se
# difference-adjusted intervals
part1$low85 <- part1$est + qnorm((1-min(vals1))/2)*part1$se
part1$up85 <- part1$est + qnorm(1-(1-min(vals1))/2)*part1$se
part1$low58 <- part1$est + qnorm((1-min(vals2))/2)*part1$se
part1$up58 <- part1$est + qnorm(1-(1-min(vals2))/2)*part1$se
# return value of function
part1
}
 
# obtain data
fac.vals <- make.data5(.75,.58,lme.model2)
 
# rewrite group panels function so that it draws two difference-adjusted intervals
 
 my.panel2 <- function(x, y, subscripts, col, pch, group.number, ...) {
#subset variables for the current panel
low95 <- fac.vals$low95[subscripts]
up95 <- fac.vals$up95[subscripts]
low85 <- fac.vals$low85[subscripts]
up85 <- fac.vals$up85[subscripts]
low58 <- fac.vals$low58[subscripts]
up58 <- fac.vals$up58[subscripts]
myjitter <- c(-.05,.05)
col2 <- c('tomato','grey60')
#95% confidence interval
panel.arrows(x+myjitter[group.number], low95, x+myjitter[group.number], up95, angle=90, code=3, length=.05, col=col)
#difference-adjusted confidence interval
panel.segments(x+myjitter[group.number], low85, x+myjitter[group.number], up85, col=col2[group.number], lineend=1, lwd=5)
panel.segments(x+myjitter[group.number], low58, x+myjitter[group.number], up58, col='dodgerblue1', lineend=1, lwd=8)
# connect means with line segments
panel.lines(x+myjitter[group.number], y, col=col)
# add point estimates for the means
panel.xyplot(x+myjitter[group.number], y, col='white', pch=16)
panel.xyplot(x+myjitter[group.number], y, col=col, pch=pch, ...)
# grid lines
panel.abline(h=seq(0,70,2), col='lightgrey', lty=3)
}
 
# generate graph
xyplot(est~L|factor(P, level=c('absent','present'), 
labels=paste('P = ', c('absent','present'), sep='')), 
xlab ='L', ylab='Estimated mean', groups=N, data=fac.vals, 
col=c(2,1), pch=c(1,16), prepanel=prepanel.ci2, ly=fac.vals$low95, 
uy=fac.vals$up95, layout=c(2,1), panel.groups=my.panel2, panel="panel.superpose",
key=list(x=.05, y=.82, corner=c(0,0), points=list(pch=c(16,1), col=1:2, cex=.9), 
text=list(c('N = 1','N = 0'), cex=.8)))
 
# we do not need both the .58 and the .75 intervals
# use just .58 interval
fac.vals <- make.data4(.58,lme.model2)
 
xyplot(est~L|factor(P, level=c('absent','present'), 
labels=paste('P = ', c('absent','present'), sep='')), 
xlab ='L', ylab='Estimated mean', groups=N, data=fac.vals, 
col=c(2,1), pch=c(1,16), prepanel=prepanel.ci2, ly=fac.vals$low95, 
uy=fac.vals$up95, layout=c(2,1), panel.groups=my.panel, panel="panel.superpose",
key=list(x=.05, y=.82, corner=c(0,0), points=list(pch=c(16,1), col=1:2, cex=.9), 
text=list(c('N = 1','N = 0'), cex=.8)))
 
#### repeated measures design ####
 
turtles <- read.table('ecol 563/turtles.txt', header=T)
turtles
 
# reorder the treatment variable
turtles$treat <- factor(turtles$treat, levels=c('fed','fast10','fast20'))
levels(turtles$treat)
 
# separate groups panel graph
xyplot(plasma~treat|sex, groups=subject, data=turtles, type='o')
 
out2 <- lme(plasma~treat*sex, random=~1|subject, data=turtles)
anova(out2)
out1 <- lme(plasma~treat+sex, random=~1|subject, data=turtles)
anova(out1)
summary(out1)
# we saw last time using MCMC that sex should be retained
 
 
# panel function version of separate groups plot
xyplot(plasma~treat|sex, groups=subject, data=turtles, type='o', panel=function(x, y, groups,...) {
panel.superpose(x,y,groups,...)})
 
fems <- predict(out1, newdata=data.frame(treat=levels(turtles$treat), sex=rep('female',3)), level=0)
fems
males <- predict(out1, newdata=data.frame(treat=levels(turtles$treat), sex=rep('male',3)), level=0)
males
 
data.frame(treat=levels(turtles$treat), sex=rep('male',3))
 
# add marginal means separately for each sex
xyplot(plasma~treat|sex, groups=subject, data=turtles, type='o', 
panel=function(x, y, subscripts, groups,...) {
 panel.superpose(x, y, subscripts, groups,...)
 sex <- turtles$sex[subscripts][1]
 if(sex=='male') panel.lines(1:3, males, type='o', pch=8, lwd=2, col=1) else
panel.lines(1:3, fems, type='o', pch=8, lwd=2, col=1)
})
 
# change individual trajectory colors to gray70
xyplot(plasma~treat|sex, groups=subject, data=turtles, type='o', 
panel=function(x, y, subscripts, groups,...) {
 panel.superpose(x, y, subscripts, groups, col='grey70',...)
 sex <- turtles$sex[subscripts][1]
 if(sex=='male') panel.lines(1:3, males, type='o', pch=8, lwd=2, col=1) else
panel.lines(1:3, fems, type='o', pch=8, lwd=2, col=1)
})